home *** CD-ROM | disk | FTP | other *** search
/ Languguage OS 2 / Languguage OS II Version 10-94 (Knowledge Media)(1994).ISO / language / parallax / more_exa.tar / more / seidel.p < prev    next >
Text File  |  1993-10-13  |  2KB  |  97 lines
  1. SYSTEM gauss_seidel;
  2. (* modified gauss-seidel iterative method *)
  3. (* see: Akl or Leighton                   *)
  4. CONST n        = 3;
  5.       eps      =   0.01;
  6.       max_iter = 100;
  7.       init_val =   1.0;
  8.  
  9. TYPE  matrix   = ARRAY [1..n],[1..n+1] OF REAL;
  10.       line     = ARRAY [1..n] OF REAL;
  11.  
  12. CONFIGURATION grid[0..n-1],[0..n];
  13. CONNECTION diag : grid[i,j] -> grid[i,i].diag;
  14.            col  : grid[i,i] -> grid[0..n-1,i].col;
  15.  
  16. SCALAR a      : matrix;
  17.        i,j    : INTEGER;
  18.        aborted: BOOLEAN;
  19.  
  20.  
  21. PROCEDURE input (SCALAR VAR a : matrix);
  22. (* reading input for A and b from terminal *)
  23. SCALAR i,j   : INTEGER;
  24. BEGIN
  25.   WriteString("Enter value of matrix A and vector b line by line :");
  26.   WriteLn;
  27.   FOR i := 1 TO n DO
  28.     FOR j := 1 TO n+1 DO ReadReal(a[i,j]) END;
  29.   END;
  30. END input;
  31.  
  32. PROCEDURE print_res(SCALAR regular: BOOLEAN; SCALAR res: line);
  33. SCALAR i: INTEGER;
  34. BEGIN
  35.   IF NOT regular THEN
  36.     WriteString("iteration aborted after max loops process ");
  37.     WriteLn;
  38.   ELSE
  39.     WriteString("Results:"); WriteLn;
  40.     FOR i:=1 TO n DO
  41.       WriteFixPt(res[i],20,7);
  42.       WriteLn;
  43.     END;
  44.   END;
  45. END print_res;
  46.  
  47.  
  48. PROCEDURE solve (SCALAR sa: matrix);
  49. SCALAR iter           : INTEGER;
  50.        regular        : BOOLEAN;
  51.        res            : line;
  52.        diff           : REAL;
  53. VECTOR a,b, x,old,nsum: REAL;
  54.   
  55. BEGIN
  56.   LOAD(a,sa);
  57.   regular := FALSE;
  58.   PARALLEL
  59.     x   := init_val;
  60.     old := init_val;
  61.     IF DIM2=n THEN SEND grid.diag(a) TO grid.diag(b) END;
  62.     iter := 0;
  63.     REPEAT
  64.       IF DIM1 <> DIM2 THEN 
  65.         IF DIM2<n THEN SEND grid.diag(a*x) TO grid.diag(nsum) REDUCE.SUM END
  66.       ELSE (* diagonal *)
  67.         x := 0.5 * x + 0.5 * (b - nsum) / a;
  68.         (* x := (b - nsum) / a; *)
  69.         diff := REDUCE.SUM( ABS(x-old) );
  70.         IF diff < eps THEN
  71.           regular := TRUE;
  72.           STORE(x,res);
  73.         ELSE
  74.           old := x;
  75.           SEND grid.col(x) TO grid.col(x);
  76.           INC(iter);
  77.         END; (* if diff *)
  78.       END; (* if DIM *)
  79.     UNTIL regular OR (iter = max_iter);
  80.   ENDPARALLEL;
  81.   print_res(regular, res);
  82. END solve;
  83.  
  84.  
  85. BEGIN (* main *)
  86.   input(a);
  87.   aborted := FALSE;
  88.   FOR i:=1 TO n DO  (* check *)
  89.     IF a[i,i] = 0.0 THEN 
  90.       WriteString("aborted: a[i,i] = 0.0");
  91.       WriteLn;
  92.       aborted := TRUE;
  93.     END
  94.   END;
  95.   IF NOT aborted THEN solve(a) END;
  96. END gauss_seidel.
  97.